Distinct pulmonary and systemic effects of dexamethasone in severe COVID-19

Summary Dexamethasone is the standard of care for critically ill patients with COVID-19, but the mechanisms by which it decreases mortality and its immunological effects in this setting are not understood. We performed bulk and single-cell RNA sequencing of the lower respiratory tract and blood, and plasma cytokine profiling to study the effect of dexamethasone on systemic and pulmonary immune cells. We find decreased signatures of antigen presentation, T cell recruitment, and viral injury in patients treated with dexamethasone. We identify compartment- and cell- specific differences in the effect of dexamethasone in patients with severe COVID-19 that are reproducible in publicly available datasets. Our results highlight the importance of studying compartmentalized inflammation in critically ill patients.


Summary 23
Dexamethasone is the standard of care for critically ill patients with COVID- 19,24 but the mechanisms by which it decreases mortality and its immunological effects in this 25 setting are not understood. We performed bulk and single-cell RNA sequencing of the 26 lower respiratory tract and blood, and plasma cytokine profiling to study the effect of 27 dexamethasone on systemic and pulmonary immune cells. We find decreased 28 signatures of antigen presentation, T cell recruitment, and viral injury in patients treated 29 with dexamethasone. We identify compartment-and cell-specific differences in the 30 effect of dexamethasone in patients with severe COVID-19 that are reproducible in 31 publicly available datasets. Our results highlight the importance of studying 32 compartmentalized inflammation in critically ill patients. Here, we use single-cell RNA sequencing to study peripheral blood and tracheal 57 aspirate (TA) from a multi-center observational cohort of patients with COVID-19 before 58 and after dexamethasone became standard of care, using data generated as part of the 59 COMET and IMPACC studies. 9, 10 We integrate this data with cytokine and gene 60 expression data from blood and compare it to two publicly available datasets. We 61 identify several cell-specific differences in the pulmonary and systemic effects of 62 dexamethasone in mechanically ventilated patients with COVID-19 ARDS, many of 63 which were reproducible in the external datasets. Through receptor-ligand analysis we 64  In order to compare systemic and tissue-specific effects of dexamethasone 187 treatment, we examined single-cell RNA sequencing data from both whole blood and 188 TA from patients treated with or without dexamethasone. We evaluated whole blood 189 (WB) scRNA-seq data from 7 Dex and 3 NoDex, and TA scRNA-seq data from 10 Dex 190 and 7 NoDex patients ( Figure 3A, 3B). A single data processing pipeline was used to 191 align, harmonize, and cluster data and identify cell types from both compartments 192 ( Figure 3C, 3D), as well as examine the cell-specific effect of dexamethasone ( Figure  193   3E, 3F). Notably, while we include in our gene expression and pathway analysis the 194 cells that are identified as neutrophils, we excluded them from our comparisons of cell 195 type abundance because their proportions were highly discordant with complete blood 196 count results of absolute neutrophil count per white blood cell count (Extended Data 197   Table 3; 202 Supplementary File 1). The greatest concordance across compartments appeared in 203 neutrophil differential gene expression (R = 0.5; Figure 3G). Dex subjects exhibited 204 decreases in expression of the S100A family of proinflammatory genes in neutrophils in 205 both lungs and blood. In contrast, gene expression in T cell subsets was highly 206 discordant across compartments (Tregs R = 0.03; CD4 T cells R = 0.05, CD8 T cells R 207 = -0.01). The greatest shared significant difference across anatomical sites in CD4 and 208 CD8 T cells was in the expression of FKBP5 (log 2 fold-difference 0.49 and 0.39, and 209 adj. p-value 0.023 and 0.058 for CD4 and CD8 T cells, respectively), which is a 210 canonical transcriptomic marker of glucocorticoid receptor activity. 18 211 In order to assess consistency and reproducibility of our analysis, we also 212 analyzed two external single-cell RNA-seq datasets using this same pipeline: Sinha et 213 al similarly generated scRNA-seq on whole blood to examine the role of neutrophils in 214 COVID-19 and responsiveness to dexamethasone in an observational cohort of 13 215 patients (5 Dex/ 8 NoDex) 7 ; and Liao et al acquired bronchoalveolar lavage (BAL) 216 samples from 6 COVID-19 patients 19 , a subset of whom were treated with the 217 corticosteroid methylprednisolone (4 methylprednisolone, 2 no-methylprednisolone). 218 Immune cell composition was similar per compartment in external datasets (Extended 219 Data Figure 5). 220 To assess whether the effects of dexamethasone were reproducible across 221 datasets, we performed tested for enrichment of pathways in the Reactome dataset that 222 were detected across blood datasets ( Figure 4A, Extended Data Figure 6) and lung   In contrast, when examining our lung datasets, we observed reproducible but 262 often discordant effects with what was observed in blood, most strikingly an elevation in 263 interferon signaling and response in influenza-related genes in T cell subsets and NK 264 cells in Dex patients that was not observed (interferon) or decreased (influenza) in the 265 blood single-cell datasets ( Figure 4B). Interferon signaling was, as expected, lower in 266 healthy controls than in COVID-19 patients (column 3). Discordant effects also included 267 pathways related to translation and cellular responses to starvation in CD4 T cells, 268 which appeared higher in lung but lower in blood in Dex patients. Concordant effects 269 across compartments were not detectable. given prior studies associating SELPLG with murine lung injury and higher risk for non-287 COVID-19 ARDS in humans. Similar effects were also observed in blood, but the effect 288 was much smaller in magnitude than in TA samples and statistically insignificant. 289 Dexamethasone was associated with additional differences in whole blood that suggesting dexamethasone may be contributing to a restoration toward a healthy 294 phenotype ( Figure 5D). The collagen and annexin pathways were more active in 295 NoDex subjects, and activity of these pathways in Dex subjects was comparable to interferon-gamma concentrations were also lower in patients treated with 340 dexamethasone. In contrast to IP-10 and IFN-gamma protein levels, interferon-341 stimulated genes were markedly upregulated in dexamethasone-treated patients in our 342 integrative analysis. The discordance between interferon levels from protein biomarker 343 data and the enrichment of interferon-related genes may reflect steroid-resistant ISG 344 pathways remaining active in these patients, which may explain the efficacy of 345 JAK/STAT inhibition in patients treated with steroids 28 . We also found higher levels of 346 Ang-1, and lower concentrations of its antagonist, Ang-2, were associated with dexamethasone treatment. An increased ratio of Ang-2 to Ang-1 reflects endothelial 348 injury 29 , and is associated with mortality in patients with ARDS due to COVID-19 and This study has several strengths. We selected subjects from a deeply 413 phenotyped observational cohort and integrated multiple assays to identify 414 compartment-and cell-specific differences in the responses to dexamethasone. We 415 build on prior studies by examining both the systemic and pulmonary biology of COVID-416 19 together, which provides more complete insight into the pathophysiology of critical 417 respiratory illness. We used mixed effects modeling to compare single cell RNA 418 expression, which addresses the pseudo-replication bias present in prior clinical single 419 cell studies and produces more conservative and reproducible estimates of differential 420 gene expression. Our findings extend our understanding of corticosteroids in critical 421 respiratory illnesses, at the gene, protein and cellular levels. Future studies using similar 422 methods can assess whether these observations are generalizable to patients with 423 other critical illness syndromes, such as sepsis or ARDS. 424 This study also has some limitations. COMET is an observational study, so 425 treatment with dexamethasone was not randomly assigned, and we cannot rule out 426 confounding by other unobserved variables that also changed during the study period.  The soluble plasma cytokines were quantified using the Luminex multiplex platform 544 (Luminex, Austin TX) as described previously. 10 Briefly, the analytes were quantified 545 using the Luminex multiplex platform with custom-developed reagents (R&D Systems, 546 Minneapolis, MN), as described in detail 43 or single-plex ELISA (R&D Systems, 547 Minneapolis, MN). The quantified analytes were read on MAGPIX® instrument and the 548 raw data was analyzed using the xPONENT® software. Analytes quantified using 549 single-plex ELISA were read using optical density. Values outside the lower limit of 550 detection were imputed using 1/3 of the lower limit of the standard curve for analytes 551 quantified by Luminex and 1/2 of the lower limit of the standard curve for analytes 552 quantified by ELISA. 553

Bulk RNA sequencing of PBMCs 555
The bulk RNA sequencing library preparation for PBMC was performed using SMART-556 Seq Low Input protocol as described here. 44 Briefly, RNA was extracted from 2.5 x 105 557 PBMCs using the Quick-RNA MagBead Kit (Zymo) with DNase digestion. RNA quality 558 was assessed using a Fragment Analyzer (Agilent) and 10ng RNA was used to 559 synthesize full length cDNA using the SMART-Seq v4 Ultra Low Input RNA Kit (Takara 560 Bio). The cDNA was purified using bead cleanup, followed by library preparation using 561 Nextera XT kit (Illumina). Libraries were validated on a Fragment Analyzer (Agilent), 562 pooled at equimolar concentrations, and sequenced on an Illumina NovaSeq6000 563 (Emory) at 100 bp paired-end read length targeting ~25 million reads per sample. 564 565

Single-cell RNA sequencing of TA and WB 566
The single cell RNA sequencing of TA and WB samples was performed as described 567 previously. 10,45 Briefly, the TA samples were transported to a BSL-3 laboratory, 3 mL of 568 TA was dissociated using 50 µg/mL collagenase type 4 (Worthington), and 0.56 ku/mL 569 of Dnase I (Worthington). The single-cells were collected by centrifugation and counted, 570 and the CD45-positive cells were enriched using MojoSort Human CD45 beads 571 (Biolgenend) and counted again before library preparation. The scRNA-seq of whole 572 blood was performed to preserve granulocytes. Briefly, the peripheral blood was 573 Gene counts were generated using the nf-core rnaseq pipeline v3.3 (https://nf-595 co.re/rnaseq) and Salmon-generated counts were used for the analyses. 596 For the analysis of bulk gene expression data, the R package DESeq2 (v1.28.1) was 597 used. Age and sex were included as covariates in the model. The log fold-change 598 values were shrunk using the apeglm algorithm. A 0.1 threshold on adjusted p-values 599 was used to identify differentially expressed genes. Gene set enrichment analysis was 600 performed with the fgsea package (v.14.0) and the REACTOME gene set database. The BCL files from sequencer were demultiplexed into individual libraries using 619 mkfastqs command in Cellranger 3.0.1 suite of tools (https://support.10xgenomics.com). 620 The feature-barcode matrices were obtained for each library by aligning the WB raw 621 FASTQ files to GRCh38 reference genome (annotated with Ensembl v85) and TA raw 622 FASTQ files to GRCh38 + SARS-CoV-2 reference genome using Cellranger count. The 623 raw feature-barcode matrices were loaded into Seurat 4.0.3, and cell barcodes with minimum of 100 features were retained in order to remove the droplets lacking cells. 625 The features that were detected in less than 3 barcodes were removed. Our dataset 626 contained three samples that were multiplexed for 10X library preparation and the rest 627 were processed individually. For the samples that were processed individually, the 628 heterotypic doublets were detected using DoubletFinder 46 by matching each cell with 629 artificially synthesized doublets. We used 35 PCs, pN=0.25 and sct=TRUE in 630 DoubletFinder. An optimal pK value (PC neighborhood size used to compute pANN) 631 was determined for each sample separately using find.pK function as suggested by the 632 authors. We approximated the doublet rate as 7% based on 10X's recommendation for 633 the expected doublets when 15,000 cells were loaded on the 10X handler 634 (https://kb.10xgenomics.com/hc/en-us/articles/360001378811). DoubletFinder requires 635 cell annotations to determine the rate of heterotypic doublets. We clustered the cell 636 barcodes using Louvain clustering and the cluster labels were used as cell annotations. 637 We removed the heterotypic doublets and subjected the remaining barcodes for further 638 quality control. 639 Our dataset contained three samples that were multiplexed, for which the filtered count 640 data for singlets were obtained from GSE163668. 10 The authors used Demuxlet 47 to 641 demultiplex the samples and to identify inter-sample doublets, and DoubletFinder to 642 identify heterotypic doublets. Single cells with greater than 50,000 unique RNA 643 molecules, fewer than 150 or greater than 8000 features, greater than 15% 644 mitochondrial content or greater than 60% ribosomal content were removed. The cell 645 cycle state of each cell was assessed using a published set of genes associated with 646 various stages of human mitosis. 48 647 The WB data from healthy controls was obtained from GSE163668, 10 the external 648 validation WB data from COVID-19 patients from GSE157789 7 and the external 649 validation bronchoalveolar lavage (BAL) fluid data from GSE145926. 19 The same data 650 processing strategy was used for these datasets as for our datasets described above. 651 652

Data integration and UMAP generation: 653
There was a substantial heterogeneity between samples within treatment groups, most 654 likely due to technical variations introduced during the library preparation that spanned 655 over months. Even if this heterogeneity is due to biological differences, this 656 heterogeneity could cause substantial issues in mapping same cell types across 657 samples. To account for this, we integrated the samples using Seurat's CCA integration 658 approach (FindIntegrationAnchors and IntegrateData functions), 49 while treating each 659 sample as its own batch. The integrated data was scaled while regressing out feature 660 counts, RNA counts, mitochondrial percentage, ribosomal percentage and cell states. 661 After reducing the data to lower dimensions (PCs), 30 PCs were used for UMAP 662 generation. The CCA integrated data was used only for generating UMAPs. All follow-663 up analyses were performed using the non-integrated data. Each tissue was processed 664 separately. 665 666 Single-cell annotation: 667 Automated cell annotation was performed using SingleR. 50 We mapped the log-668 normalized expression data against a reference expression dataset from ENCODE 669 To study the cell-type-specific effects of dexamethasone in whole blood and TA 682 samples, we compared gene expression between Dex and NoDex samples within each 683 tissue for every cell-type separately. The differential expression analysis was performed 684 objects. For each cell type, the Seurat object was subsetted to keep single-cell 688 expression data for that cell type, the subsetted object was converted to 689 SingleCellExperiment object, and the RNA raw counts were normalized for the library 690 size (i.e. divided each count by total number of UMIs per cell and multiply by the mean 691 of the number of UMIs per cell across all cells) and log 2 transformed with pseudocount 692 of 1. To remove the highly sparse data, only genes with non-zero counts in at least 5% 693 cells in at least one condition were retained. Finally, the zlm function was used to 694 identify the differentially expressed genes between Dex and NoDex samples. We 695 accounted for the number of detected genes per cell in the model. Since the numbers of 696 cells per patient are often very different, the differential analysis is often biased toward 697 the patient with the largest number of cells. To account for this bias, we used patient ids 698 as a random effect. Additionally, we used the following parameters in zlm function: 699 method='glmer', ebayes = F, strictConvergence = FALSE, fitArgsD = list(nAGQ = 0). 700 Finally, the P values were corrected for multiple testing using FDR. We performed CellChat analysis 53 to identify ligand-receptor pairs that display 712 differential interaction strength between cells from Dex, NoDex and healthy groups. The 713 Seurat objects were subsetted to include the cell types that had more than 100 cells in 714 all conditions within that tissue. Specifically, for TA data, the cell types with more than 715 100 cells in both Dex and NoDex were retained, and for blood data, the cell types with 716 more than 100 cells in all groups (Dex (COMET), NoDex (COMET), healthy (COMET), 717 Dex (Sinha et al) and NoDex (Sinha et al)) were retained. The CellChat objects were 718 first created for each group (condition) of cells separately using createCellChat() 719 function, with Seurat's normalized RNA data as input data. The over expressed genes 720 and interactions were identified based on the CellChat database of human ligand-721 receptor pairs, and the expressed data were projected on the protein-protein interaction 722 network. Finally, the communication probabilities were calculated, the communications 723 based on less than 10 cells were discarded, aggregated network were calculated by 724 summarizing the communication probability, and saved as individual RDS files for each 725 condition. Pairs of conditions, for example TA Dex and TA NoDex, were compared 726 using rankNet to rank signaling networks based on the information flow. We used this 727 information flow to find ligand-receptor pairs that exhibit significant difference in 728 predicted interaction strength between the conditions. 729

Statistics 731
The p-values were corrected for multiple testing using Benjamini-Hochberg method, 732 which controls for the false-discovery rate (FDR